/*
 * CoupledBoundaryCarbonBasedAblationMaterial.C
 *
 *  Created on: Jan 3, 2020
 *      Author: liuxiao
 */

#include "CoupledBoundaryCarbonBasedAblationMaterial.h"
#include<iostream>
#include "libmesh/quadrature.h"
#include "libmesh/system.h"
#include "libmesh/radial_basis_interpolation.h"
using std::cout;
using std::endl;
registerMooseObject("TrilobitaApp", CoupledBoundaryCarbonBasedAblationMaterial);
template<>
InputParameters validParams<CoupledBoundaryCarbonBasedAblationMaterial>()
{
	InputParameters params = validParams<Material>();
	params.addParam<Real>("M_p", 0.44, " average molar mass of pyrolysis gas,kg/mol");
	params.addParam<Real>("mp", 0, " mass flux of pyrolysis gas,kg/m**2/s");
	params.addRequiredCoupledVar("Tw", "Coupled wall temperature");
	params.addParam<int>("gamma", 0, "to choose laminar 0 or turbulence 1");
	params.addParam<int>("scala_model", 1, "to choose scala model,0=slow reaction,1=fast reaction");
	params.addParam<Real>("bw_res", 1e-4, " residul to compute bw");

	return params;
}


CoupledBoundaryCarbonBasedAblationMaterial::CoupledBoundaryCarbonBasedAblationMaterial(const InputParameters& parameters) :
	Material(parameters),
	_Pw_C1(declareProperty<Real>("Pw_C1")),
	_Pw_C2(declareProperty<Real>("Pw_C2")),
	_Pw_C3(declareProperty<Real>("Pw_C3")),
	_Pw_C4(declareProperty<Real>("Pw_C4")),
	_Pw_C5(declareProperty<Real>("Pw_C5")),
	_Pw_O2(declareProperty<Real>("Pw_O2")),
	_Pw_O(declareProperty<Real>("Pw_O")),
	_Pw_CO(declareProperty<Real>("Pw_CO")),
	_Pw_CO2(declareProperty<Real>("Pw_CO2")),
	_Pw_N2(declareProperty<Real>("Pw_N2")),
	_Pw_N(declareProperty<Real>("Pw_N")),
	_Pw_CN(declareProperty<Real>("Pw_CN")),
	_Pw_C2N(declareProperty<Real>("Pw_C2N")),
	_Cw_C1(declareProperty<Real>("Cw_C1")),
	_Cw_C2(declareProperty<Real>("Cw_C2")),
	_Cw_C3(declareProperty<Real>("Cw_C3")),
	_Cw_C4(declareProperty<Real>("Cw_C4")),
	_Cw_C5(declareProperty<Real>("Cw_C5")),
	_Cw_O2(declareProperty<Real>("Cw_O2")),
	_Cw_O(declareProperty<Real>("Cw_O")),
	_Cw_CO(declareProperty<Real>("Cw_CO")),
	_Cw_CO2(declareProperty<Real>("Cw_CO2")),
	_Cw_N2(declareProperty<Real>("Cw_N2")),
	_Cw_N(declareProperty<Real>("Cw_N")),
	_Cw_CN(declareProperty<Real>("Cw_CN")),
	_Cw_C2N(declareProperty<Real>("Cw_C2N")),
	_M_bar(declareProperty<Real>("M_bar")),
	_Bc(declareProperty<Real>("Bc")),
	_Bp(declareProperty<Real>("Bp")),
	_Bw(declareProperty<Real>("Bw")),
	_mc(declareProperty<Real>("mc")),
	_recession_rate(declareProperty<Real>("recession_rate")),
	_mw(declareProperty<Real>("mw")),
	_injection_coff(declareProperty<Real>("injection_coff")),
	_hr(getMaterialProperty<Real>("recovery_enthalpy")),
	_Pe_atm(getMaterialProperty<Real>("Pe_atm")),
	_Pe_O2(getMaterialProperty<Real>("Pe_O2")),
	_Pe_N2(getMaterialProperty<Real>("Pe_N2")),
	_M_e(getMaterialProperty<Real>("M_e")),
	_Rhos(getADMaterialProperty<Real>("solid_density")),
	_qc_0K(getMaterialProperty<Real>("heatflux_0K")),
	_M_p(getParam<Real>("M_p")),
	_mp(getParam<Real>("mp")),
	_Tw(coupledValue("Tw")),
	_gamma(getParam<int>("gamma")),
	_scala_model(getParam<int>("scala_model")),
	_bw_res(getParam<Real>("bw_res")),
	_Kp_C1(0),
	_Kp_C2(0),
	_Kp_C3(0),
	_Kp_C4(0),
	_Kp_C5(0),
	_Kp_O(0),
	_Kp_CO2(0),
	_Kp_CN(0),
	_Kp_C2N(0),
	_Kp_N(0),
	_M_C(0.012),
	_M_C1(0.012),
	_M_C2(0.024),
	_M_C3(0.036),
	_M_C4(0.048),
	_M_C5(0.060),
	_M_O2(0.032),
	_M_O(0.016),
	_M_CO(0.028),
	_M_CO2(0.044),
	_M_N2(0.028),
	_M_N(0.014),
	_M_CN(0.026),
	_M_C2N(0.046),
	_R(8.314)
{

}

void CoupledBoundaryCarbonBasedAblationMaterial::compute_Kp(Real T)
{
	if (T < 300)
	{
		T = 300;
	}
	_Kp_C1 = pow(10, (8.12 - 3.7219e4 / T));
	_Kp_C2 = pow(10, (9.64 - 4.2713e4 / T));
	_Kp_C3 = pow(10, (9.78 - 4.0481e4 / T));
	_Kp_C4 = pow(10, (10.013 - 4.889e4 / T));
	_Kp_C5 = pow(10, (10.464 - 4.927e4 / T));
	_Kp_O = pow(10, (3.518 - 13392 / T));
	_Kp_CO2 = 10 * pow(10, (-4.482 + 1.4608e4 / T));
	_Kp_CN = pow(10, (5.029 - 22288 / T));
	_Kp_C2N = pow(10, (6.572 - 28450 / T));
	_Kp_N = pow(10, (3.542 - 25176 / T));
}

Real CoupledBoundaryCarbonBasedAblationMaterial::compute_Ci_from_P(Real P, Real Mi, Real M_bar, int qp)
{
	Real C_i;
	C_i = (P / _Pe_atm[qp]) * (Mi / M_bar);
	return C_i;
}

Real CoupledBoundaryCarbonBasedAblationMaterial::compute_P_Ci(int i, Real Bw, Real phi, Real M_bar, Real T, int qp)
{
	std::vector<Real> P_Ci;
	std::vector<Real> beta;
	std::vector<Real> alpha;
	std::vector<Real> P_Ceq;
	std::vector<Real> M;
	Real Phol;
	Real Ahol;
	beta.push_back(0);
	beta.push_back(0.24);
	beta.push_back(0.5);
	beta.push_back(0.023);
	beta.push_back(0.25);
	beta.push_back(0.0019);
	alpha.push_back(0);
	alpha.push_back(beta[1] * M_bar * _Pe_atm[qp] / sqrt(2 * pi * _R * _M_C1 * T));
	alpha.push_back(beta[2] * M_bar * _Pe_atm[qp] / sqrt(2 * pi * _R * _M_C2 * T));
	alpha.push_back(beta[3] * M_bar * _Pe_atm[qp] / sqrt(2 * pi * _R * _M_C3 * T));
	alpha.push_back(beta[4] * M_bar * _Pe_atm[qp] / sqrt(2 * pi * _R * _M_C4 * T));
	alpha.push_back(beta[5] * M_bar * _Pe_atm[qp] / sqrt(2 * pi * _R * _M_C5 * T));
	P_Ceq.push_back(0);
	P_Ceq.push_back(_Kp_C1);
	P_Ceq.push_back(_Kp_C2);
	P_Ceq.push_back(_Kp_C3);
	P_Ceq.push_back(_Kp_C4);
	P_Ceq.push_back(_Kp_C5);
	M.push_back(0);
	M.push_back(_M_C1);
	M.push_back(_M_C2);
	M.push_back(_M_C3);
	M.push_back(_M_C4);
	M.push_back(_M_C5);
	Phol = phi * _qc_0K[qp] / _hr[qp];
	Ahol = M_bar * _Pe_atm[qp];
	for (unsigned int i(0); i < 6; ++i)
	{
		if (i == 0)
		{
			P_Ci.push_back(0);
		}
		else
		{
			P_Ci.push_back(P_Ceq[i] / (1 + (1 + Bw) * Phol / (alpha[i] * Ahol / M[i])));
		}

	}

	return P_Ci[i];

}

Real CoupledBoundaryCarbonBasedAblationMaterial::compute_Pw_O2(Real Pw_CO, Real Pw_CO2, Real Pw_O, Real Bw, Real phi, Real M_bar, Real T, int qp)
{
	Real A1;
	Real B1;
	Real alpha1;
	Real A1_p;
	Real B1_p;
	Real alpha1_p;
	Real Phol;
	Real a;
	Real b;
	Real c;
	Real X;
	Real P_CO_CO2;
	Real P_O_CO2;
	Real alch;
	int f;
	if (T < 300)
	{
		T = 300;
	}
	if (_scala_model == 0)
	{
		A1 = 2.18e5;
		B1 = 21558;
		alpha1 = A1 * exp(-B1 / T);
	}
	else if (_scala_model == 1)
	{
		A1 = 3.29e9;
		B1 = 22144;
		alpha1 = A1 * exp(-B1 / T);
	}
	else
	{
		mooseError("wrong scala_model value,must be 0 or 1");
	}
	A1_p = 4.88e10;
	B1_p = 42500;
	alpha1_p = A1_p * exp(-B1_p / T);
	//   if(T<=3300)
	//   {
	f = 0;
	//   }
	//   else
	//   {
	//	   f=1;
	//   }

	if (abs(Pw_CO2 < 1e-10))
	{
		Pw_CO2 = 1e-10;
	}
	P_CO_CO2 = Pw_CO / Pw_CO2;
	Pw_O = Pw_O / Pw_CO2;
	alch = (2 + P_CO_CO2 + Pw_O) * _M_O2 / (2 + 2 * P_CO_CO2) / _M_C;
	Phol = phi * _qc_0K[qp] / _hr[qp];
	a = Phol * _M_O2 * (1 + Bw) / M_bar + f * alpha1_p * alch;
	b = (1 - f) * alpha1 * alch;
	c = -Phol * _M_O2 * _Pe_O2[qp] / M_bar;
	if ((pow(b, 2) - 4 * a * c) > 0 && abs(a) > 1e-32)
	{
		X = (-b + sqrt(pow(b, 2) - 4 * a * c)) / (2 * a);
	}
	else
	{
		X = 0;
		mooseError("Pw_O2 has no solution!");
	}
	Real mc1 = (1 - f) * alpha1 * X + f * alpha1_p * pow(X, 2);

	return pow(X, 2);

}
Real CoupledBoundaryCarbonBasedAblationMaterial::compute_Pw_O(Real Pw_O2)
{
	Real Pw_O;
	Pw_O = Pw_O2 * _Kp_O;
	return Pw_O;
}

Real CoupledBoundaryCarbonBasedAblationMaterial::compute_Pw_CO(Real Pw_O2, Real Pw_O, Real M_bar, Real Bw, int qp)
{
	Real Ce_O;
	Real Pw_CO;
	Ce_O = compute_Ci_from_P(_Pe_O2[qp], _M_O2, _M_e[qp], qp);
	Pw_CO = (Ce_O * M_bar * _Pe_atm[qp] / (1 + Bw) / _M_O - Pw_O - 2 * Pw_O2) / (1 + 2 * sqrt(Pw_O2) * _Kp_CO2);
	return Pw_CO;
}

Real CoupledBoundaryCarbonBasedAblationMaterial::compute_Pw_CO2(Real Pw_CO, Real Pw_O2)
{
	Real Pw_CO2;
	Pw_CO2 = _Kp_CO2 * Pw_CO * sqrt(Pw_O2);
	return Pw_CO2;
}

Real CoupledBoundaryCarbonBasedAblationMaterial::compute_Pw_N2(Real Bw, Real phi, Real M_bar, int qp)
{
	Real Ce_N;
	Real a;
	Real b;
	Real c;
	Real Y;
	Ce_N = compute_Ci_from_P(_Pe_N2[qp], _M_N2, _M_e[qp], qp);
	a = 2;
	b = _Kp_N + _Kp_CN + _Kp_C2N;
	c = -_Pe_atm[qp] * M_bar * Ce_N / _M_N / (1 + Bw);
	if ((pow(b, 2) - 4 * a * c) > 0)
	{
		Y = (-b + sqrt(pow(b, 2) - 4 * a * c)) / (2 * a);

	}
	else
	{
		Y = 0;
		//mooseError("Pw_N2 has no solution!");
	}
	return pow(Y, 2);
}
Real CoupledBoundaryCarbonBasedAblationMaterial::compute_Pw_N(Real Pw_N2)
{
	Real Pw_N;
	Pw_N = _Kp_N * sqrt(Pw_N2);
	return Pw_N;
}

Real CoupledBoundaryCarbonBasedAblationMaterial::compute_Pw_CN(Real Pw_N2)
{
	Real Pw_CN;
	Pw_CN = _Kp_CN * sqrt(Pw_N2);
	return Pw_CN;
}
Real CoupledBoundaryCarbonBasedAblationMaterial::compute_Pw_C2N(Real Pw_N2)
{
	Real Pw_C2N;
	Pw_C2N = _Kp_C2N * sqrt(Pw_N2);
	return Pw_C2N;
}

Real CoupledBoundaryCarbonBasedAblationMaterial::compute_M_bar(Real Cw_C1, Real Cw_C2, Real Cw_C3, Real Cw_C4, Real Cw_C5, Real Cw_O2, Real Cw_O, Real Cw_CO, Real Cw_CO2, Real Cw_N2, Real Cw_N, Real Cw_CN, Real Cw_C2N)
{
	Real M_bar;
	M_bar = 1 / (Cw_C1 / _M_C1 + Cw_C2 / _M_C2 + Cw_C3 / _M_C3 + Cw_C4 / _M_C4 + Cw_C5 / _M_C5\
		+ Cw_O / _M_O + Cw_O2 / _M_O2 + Cw_CO / _M_CO + Cw_CO2 / _M_CO2 + Cw_N / _M_N + Cw_N2 / _M_N2\
		+ Cw_CN / _M_CN + Cw_C2N / _M_C2N);
	return M_bar;

}

Real CoupledBoundaryCarbonBasedAblationMaterial::compute_carbon_w(Real Cw_C1, Real Cw_C2, Real Cw_C3, Real Cw_C4, Real Cw_C5, \
	Real Cw_CO, Real Cw_CO2, Real Cw_CN, Real Cw_C2N)
{
	Real carbon_w;
	carbon_w = Cw_C1 + Cw_C2 + Cw_C3 + Cw_C4 + Cw_C5 + Cw_CO * _M_C / _M_CO + Cw_CO2 * _M_C / _M_CO2\
		+ Cw_CN * _M_C / _M_CN + Cw_C2N * _M_C / _M_C2N;
	return carbon_w;
}

Real CoupledBoundaryCarbonBasedAblationMaterial::compute_Bc(Real carbon_w, Real Bp)
{

	Real Bc;
	if (carbon_w > 0.999)
	{
		carbon_w = 0.999;
	}
	Bc = (1 + Bp) * carbon_w / (1 - carbon_w);
	return Bc;
}

Real CoupledBoundaryCarbonBasedAblationMaterial::compute_Bp(Real mp, Real phi, int qp)
{

	Real Bp;
	Bp = mp / (phi * _qc_0K[qp] / _hr[qp]);
	return Bp;
}
Real CoupledBoundaryCarbonBasedAblationMaterial::compute_mc(Real Bc, Real phi, int qp)
{
	Real mc;
	mc = Bc * phi * _qc_0K[qp] / _hr[qp];
	return mc;
}


Real CoupledBoundaryCarbonBasedAblationMaterial::compute_phi(Real Bc, Real M_bar, Real mc, Real mp, int gamma, Real T, int qp)
{
	Real kesi = (pow((M_bar / _M_C), 0.26) * mc + pow((M_bar / _M_p), 0.26) * mp) * _hr[qp] / _qc_0K[qp];
	Real phi;
	if (T < 3300)
	{
		if (kesi < 2.25)
		{
			phi = 1 - 0.724 * kesi + 0.13 * pow(kesi, 2);
		}
		else
		{
			phi = 0.04;
		}
	}
	else
	{
		phi = 1 - (1 - gamma) * 0.58 * (mc + mp) * _hr[qp] / _qc_0K[qp] - gamma * 2 * (mc + mp) * _hr[qp] / _qc_0K[qp];
	}

	return phi;

}

void CoupledBoundaryCarbonBasedAblationMaterial::iteration(Real T, int qp)
{
	Real Cw_C1_it = 0;
	Real Cw_C2_it = 0;
	Real Cw_C3_it = 0;
	Real Cw_C4_it = 0;
	Real Cw_C5_it = 0;
	Real Cw_O2_it = 0;
	Real Cw_O_it = compute_Ci_from_P(_Pe_O2[qp], _M_O2, _M_e[qp], qp);
	Real Cw_CO_it = 0;
	Real Cw_CO2_it = 0;
	Real Cw_N2_it = compute_Ci_from_P(_Pe_N2[qp], _M_N2, _M_e[qp], qp);
	Real Cw_N_it = 0;
	Real Cw_CN_it = 0;
	Real Cw_C2N_it = 0;
	Real Pw_C1_it = 0;
	Real Pw_C2_it = 0;
	Real Pw_C3_it = 0;
	Real Pw_C4_it = 0;
	Real Pw_C5_it = 0;
	Real Pw_O2_it = _Pe_O2[qp];
	Real Pw_O_it = 0;
	Real Pw_CO_it = 0;
	Real Pw_CO2_it = 0;
	Real Pw_N2_it = _Pe_N2[qp];
	Real Pw_N_it = 0;
	Real Pw_CN_it = 0;
	Real Pw_C2N_it = 0;
	Real carbon_w_it = 0;
	Real mc_it = 0;
	Real Bc_it0 = 0;
	Real Bc_it1 = 0;
	Real Bp_it = 0;
	Real delta_Bc_it = 0.0001;
	Real M_bar_it = _M_e[qp];
	Real phi_it = 1;
//	Real Res = 1e-4;
	Real Bw_it = 0;
	int init = 1;
	compute_Kp(T);
	while (init == 1 || abs(Bc_it1 - Bc_it0) > _bw_res)
	{  
		init = init + 1;
		if (init > 5000)
		{
		mooseError("bw iteration>5000,cannot converge");
		}
		delta_Bc_it = (Bc_it1 - Bc_it0) / 2;
		Bc_it0 = Bc_it0 + delta_Bc_it;
	//	Bc_it0 = Bc_it1;
		Bp_it = compute_Bp(_mp, phi_it, qp);
		Bw_it = Bc_it0 + Bp_it;
		Pw_O2_it = compute_Pw_O2(Pw_CO_it, Pw_CO2_it, Pw_O_it, Bw_it, phi_it, M_bar_it, T, qp);
		Pw_O_it = compute_Pw_O(Pw_O2_it);
		Pw_CO_it = compute_Pw_CO(Pw_O2_it, Pw_O_it, M_bar_it, Bw_it, qp);
		Pw_CO2_it = compute_Pw_CO2(Pw_CO_it, Pw_O2_it);
		Cw_O2_it = compute_Ci_from_P(Pw_O2_it, _M_O2, M_bar_it, qp);
		Cw_O_it = compute_Ci_from_P(Pw_O_it, _M_O, M_bar_it, qp);
		Cw_CO_it = compute_Ci_from_P(Pw_CO_it, _M_CO, M_bar_it, qp);
		Cw_CO2_it = compute_Ci_from_P(Pw_CO2_it, _M_CO2, M_bar_it, qp);
		if (T >= 3000)
		{
			Pw_C1_it = compute_P_Ci(1, Bw_it, phi_it, M_bar_it, T, qp);
			Pw_C2_it = compute_P_Ci(2, Bw_it, phi_it, M_bar_it, T, qp);
			Pw_C3_it = compute_P_Ci(3, Bw_it, phi_it, M_bar_it, T, qp);
			Pw_C4_it = compute_P_Ci(4, Bw_it, phi_it, M_bar_it, T, qp);
			Pw_C5_it = compute_P_Ci(5, Bw_it, phi_it, M_bar_it, T, qp);
			Pw_N2_it = compute_Pw_N2(Bw_it, phi_it, M_bar_it, qp);
			Pw_N_it = compute_Pw_N(Pw_N2_it);
			Pw_CN_it = compute_Pw_CN(Pw_N2_it);
			Pw_C2N_it = compute_Pw_C2N(Pw_N2_it);
			Cw_C1_it = compute_Ci_from_P(Pw_C1_it, _M_C1, M_bar_it, qp);
			Cw_C2_it = compute_Ci_from_P(Pw_C2_it, _M_C2, M_bar_it, qp);
			Cw_C3_it = compute_Ci_from_P(Pw_C3_it, _M_C3, M_bar_it, qp);
			Cw_C4_it = compute_Ci_from_P(Pw_C4_it, _M_C4, M_bar_it, qp);
			Cw_C5_it = compute_Ci_from_P(Pw_C5_it, _M_C5, M_bar_it, qp);
			Cw_N2_it = compute_Ci_from_P(Pw_N2_it, _M_N2, M_bar_it, qp);
			Cw_N_it = compute_Ci_from_P(Pw_N_it, _M_N, M_bar_it, qp);
			Cw_CN_it = compute_Ci_from_P(Pw_CN_it, _M_CN, M_bar_it, qp);
			Cw_C2N_it = compute_Ci_from_P(Pw_C2N_it, _M_C2N, M_bar_it, qp);
		}
		carbon_w_it = compute_carbon_w(Cw_C1_it, Cw_C2_it, Cw_C3_it, Cw_C4_it, Cw_C5_it, Cw_CO_it, Cw_CO2_it, Cw_CN_it, Cw_C2N_it);
		Bc_it1 = compute_Bc(carbon_w_it, Bp_it);
		mc_it = compute_mc(Bc_it1, phi_it, qp);
		M_bar_it = compute_M_bar(Cw_C1_it, Cw_C2_it, Cw_C3_it, Cw_C4_it, Cw_C5_it, Cw_O2_it, Cw_O_it, Cw_CO_it, Cw_CO2_it, Cw_N2_it, Cw_N_it, Cw_CN_it, Cw_C2N_it);
		phi_it = compute_phi(Bc_it1, M_bar_it, mc_it, _mp, _gamma, T, qp);
		Real P_CO_CO2 = Pw_CO_it / (Pw_CO2_it + 1e-10);
	}
	_Cw_C1[qp] = Cw_C1_it;
	_Cw_C1[qp] = Cw_C1_it;
	_Cw_C2[qp] = Cw_C2_it;
	_Cw_C3[qp] = Cw_C3_it;
	_Cw_C4[qp] = Cw_C4_it;
	_Cw_C5[qp] = Cw_C5_it;
	_Cw_O2[qp] = Cw_O2_it;
	_Cw_O[qp] = Cw_O_it;
	_Cw_CO[qp] = Cw_CO_it;
	_Cw_CO2[qp] = Cw_CO2_it;
	_Cw_N2[qp] = Cw_N2_it;
	_Cw_N[qp] = Cw_N_it;
	_Cw_CN[qp] = Cw_CN_it;
	_Cw_C2N[qp] = Cw_C2N_it;
	_Bc[qp] = Bc_it0;
	_Bw[qp] = _Bc[qp] + _Bp[qp];
	_M_bar[qp] = M_bar_it;
	_injection_coff[qp] = phi_it;
	_mc[qp] = _Bc[qp] * _injection_coff[qp] * _qc_0K[qp] / _hr[qp];
	_recession_rate[qp] = _mc[qp] / _Rhos[qp].value();
}


void CoupledBoundaryCarbonBasedAblationMaterial::computeProperties()
{

	for (unsigned int qp(0); qp < _qrule->n_points(); ++qp)
	{

		iteration(_Tw[qp], qp);

	}

}
